rm(list=ls())
source('param_general.R')


source('./SwitchCosts/stata_results_switch.R')
## Paper Table 2, Age Restriction vs. Honore 
DISTANCE = c(coefDistList$main,coefDistList$naive, coefDistList$naive)

SWITCHING.COST =  c(coefSwitchList$main,coefSwitchList$naive, 0)


dataPat <- readRDS(paste0(datDir,"ld_women_process_all.rds"))
dataPat = copy(dataPat[,list(maskssn,hosp_id,pat_zip,age,yearDat)])
setnames(dataPat,"yearDat","year")
setkey(dataPat, maskssn,year)
dataPat[,':='(mult_admit=.N,birthNum = .I - min(.I) + 1L),by=list(maskssn)]
dataPat[mult_admit>1,':='(hosp_id_prev=c(NA,head(hosp_id,-1)),
                      pat_zip_prev=c(NA,head(pat_zip,-1)),
                      year_prev=c(NA,head(year,-1))),
    by=list(maskssn)]
dataPat[,years_since_birth:= year - year_prev ]

PROB_NEW_MOTHER = length(dataPat[year == 2014 & birthNum ==1]$birthNum)/length(dataPat[year == 2014]$birthNum)

PROB_OLD_MOTHER = table(dataPat[year == 2014 & birthNum >1 & years_since_birth > 0]$years_since_birth)/length(dataPat[year == 2014 & birthNum > 1 & years_since_birth > 0]$birthNum)

source("./SwitchCosts/create_delta_entry.R")

dataPatPSA = copy(dataPat[pat_zip %in% zipPSA,])

## Now create cohorts ...

dataCohortPSA = copy(dataPatPSA[birthNum == mult_admit,list(maskssn,hosp_id,pat_zip,age,year,birthNum)])
setnames(dataCohortPSA,"hosp_id","hosp_id_alt")
dataCohortPSA[,hosp_id_alt := as.character(hosp_id_alt) ]
dataCohortPSA = merge(dataCohortPSA,hospitalCrosswalk,by = c("hosp_id_alt"))
dataCohortPSA[!(hospital %in% hospitalsInMarket), hosp_id_alt :="0"]
dataCohortPSA[, hospital := NULL]

setnames(dataCohortPSA,c("hosp_id_alt","year","age"),c("hosp_id_prev","year_prev","age_prev"))

## Do 2015
## Now cohort demand based on 2013 actuals
## Adjusted to future as in crete_delta_conlaw.R
demandTotal = dataPatPSA[year == 2013,.N,by=c("pat_zip")]
setnames(demandTotal,"N","nZipNew")
demandTotal[pat_zip %in% ZIP_CLOSER,nZipNew := nZipNew * POP_ADJ_CLOSER]
demandTotal[pat_zip %in% ZIP_CONSTANT,nZipNew := nZipNew * POP_ADJ_CONSTANT]
demandTotal[pat_zip %in% ZIP_FARTHER,nZipNew := nZipNew * POP_ADJ_FARTHER]
demandTotal[,newMother := round(PROB_NEW_MOTHER*round(nZipNew))]

## New Mothers
newMotherCohort = NULL
for(xi in 1: length(demandTotal$pat_zip)) {
  temp = as.data.table(list(admission= 1:demandTotal[xi,newMother],pat_zip = demandTotal[xi,pat_zip]))
  temp[,maskssn := paste0(as.character(pat_zip),"-",as.character(admission))]
  newMotherCohort = rbind(newMotherCohort,temp[,list(maskssn,pat_zip)])
}
## Old Mothers
prevMotherProb =   as.data.frame(list(years_since_birth = 1:8,prob =  unique(PROB_OLD_MOTHER)))

demandOld = as.data.table(merge(as.data.frame(demandTotal),prevMotherProb,by= NULL))
demandOld[,demand := round(prob *(1-PROB_NEW_MOTHER) * round(nZipNew))]


dataPSAMktList = list(all = copy(dataPSAMkt), naive = copy(dataPSAMkt), noSC = copy(dataPSAMkt))


for(xi in 1:3) {
  dataPSAMktList[[xi]][,xiHat := deltaHat - DISTANCE[xi] * (time_current)]
  dataPSAMktList[[xi]][hosp_id_alt != "999999",meanXi := mean(xiHat), by = c("pat_zip")]
  dataPSAMktList[[xi]][hosp_id_alt != "999999",sdXi := sd(xiHat)]
  dataPSAMktList[[xi]][,meanXi := mean(xiHat, na.rm = TRUE), by = c("pat_zip")]
  dataPSAMktList[[xi]][,sdXi := sd(xiHat, na.rm = TRUE), by = c("pat_zip")]
  
  
  dataPSAMktList[[xi]] = dataPSAMktList[[xi]][,list(hosp_id_alt,pat_zip,xiHat,time_current_new,deltaHat, meanXi, sdXi)]
  ## Add Outside Option
  dataPSAMktList[[xi]] = rbind(dataPSAMktList[[xi]],as.data.table(list(hosp_id_alt = "0",xiHat = 0, time_current_new = 0, meanXi = NA, sdXi = NA, deltaHat = 0, pat_zip = dataPSAMkt[,unique(pat_zip)])))
}


createInitialCohort = function(dataCohortInitial,dataCohortNew,demandOld,YEAR.START,YEAR.END) {
  allBirthCohorts = NULL
  
  ## Initialize Old Cohort
  cohort = copy(dataCohortInitial)
  cohort[,years_since_birth := 2015 - year_prev]
  
  for(YEAR in YEAR.START:YEAR.END) {
    cohort = merge(cohort,demandOld[,list(pat_zip,years_since_birth,demand)], by= c("pat_zip","years_since_birth"))

    pregnant_women = cohort[,sample(maskssn,demand,FALSE),by = c("pat_zip","years_since_birth")]$V1
    
    cohort[,birth_this_year := (maskssn %in% pregnant_women)]
    ## Next Cohort
    nextCohortNew = copy(dataCohortNew)[,':='(year_prev = YEAR,birthNum = 1L, years_since_birth = 1,hosp_id_prev = -1, birth_this_year = TRUE, maskssn = paste0(YEAR,"-",maskssn))]
    nextCohortOld = copy(cohort)
    nextCohortOld[,":="(demand = NULL,age_prev = NULL)]
    nextCohortOld[birth_this_year == TRUE,birthNum := birthNum + 1L]
    birthCohort = rbind(copy(nextCohortOld)[birth_this_year == TRUE,],copy(nextCohortNew)[,':='(year_prev = NA, years_since_birth = 0L)])
    birthCohort[,year := YEAR]
    nextCohortOld[birth_this_year == TRUE, ':='(year_prev = YEAR,years_since_birth = 1L)]
    nextCohortOld[birth_this_year == FALSE, ':='(years_since_birth = years_since_birth + 1L)]
    
    nextCohortOld = nextCohortOld[years_since_birth <= 8,]
    nextCohort = rbind(nextCohortNew,nextCohortOld)
    allBirthCohorts = rbind(allBirthCohorts,birthCohort)
    cohort = copy(nextCohort)
  }
  return(allBirthCohorts)
  
}

doSimulation = function(ITER,data,YEAR.START,YEAR.END,DISTANCE, SWITCHING.COST) {
  resultsList= list(mean = list(),plusHalf= list(), plusOne= list(), minusHalf = list(), minusOne = list(), old = list())
    for(xj in 1:6) {
      dataSimul = copy(data)
        for(YEAR in YEAR.START:YEAR.END) {
          dataSimul[,xiHatNew := xiHat]
          if(xj == 1) {
            dataSimul[hosp_id_alt == "999999", xiHatNew := meanXi]
          }
          if(xj == 2) {
            dataSimul[hosp_id_alt == "999999", xiHatNew := meanXi + .5 * sdXi]
          }
          if(xj == 3) {
            dataSimul[hosp_id_alt == "999999", xiHatNew := meanXi + 1 * sdXi]
          }
          if(xj == 4) {
            dataSimul[hosp_id_alt == "999999", xiHatNew := meanXi - .5 * sdXi]
          }
          if(xj == 5) {
            dataSimul[hosp_id_alt == "999999", xiHatNew := meanXi - 1 * sdXi]
          }
          if(xj == 6) {
            dataSimul = dataSimul[hosp_id_alt != "999999", ]

          }
          dataSimul[year == YEAR, detUtil := xiHatNew + DISTANCE * (time_current_new) + SWITCHING.COST * (chosenPrev == TRUE)]
          ## Outside option has no utility, no switching benefit
          dataSimul[year == YEAR & hosp_id_alt == 0, detUtil := 0]
          
          dataSimul[year == YEAR, hospShare := exp(detUtil)/sum(exp(detUtil)), by = c("maskssn","birthNum")]
          dataSimul[year == YEAR, hospShareCumSum := cumsum(hospShare), by = c("maskssn","birthNum")]
          dataSimul[year == YEAR, obs_choice := (rand<= hospShareCumSum & rand > (hospShareCumSum - hospShare))]
          dataSimul[year == YEAR,utility := log(0+sum(exp(detUtil))),by = c("maskssn","birthNum")] 
          setkey(dataSimul,maskssn,hosp_id_alt,year)
 
          dataSimul[, chosenPrevNew:= c(NA,head(obs_choice,-1)), by = list(maskssn,hosp_id_alt)]
          dataSimul[!(is.na(chosenPrevNew) | hosp_id_alt == "0"), chosenPrev := chosenPrevNew]
        }
        dataWTP1 = dataSimul[hosp_id_alt == "100167", list(share = sum(obs_choice)/.N,
                                                           share_Old = sum(obs_choice *(birthNum >1))/sum(birthNum >1),
                                                           WTP = sum(log(1/(1-hospShare))),
                                                           WTP_Old = sum(log(1/(1-hospShare))*(birthNum >1)), N = .N,
                                                           NOld = sum(birthNum >1)),
                             by = c("pat_zip","year")]
        dataWTP1[, hospital := 1]
        hosp2List = c("100038", "111527", "23960050")
        dataSimulAgg = dataSimul[hosp_id_alt %in% hosp2List, list(obs_choice = sum(obs_choice), hospShare = sum(hospShare)), by = c("maskssn","year", "birthNum", "pat_zip")]
        dataWTP2 = dataSimulAgg[, list(share = sum(obs_choice)/.N,
                                                           share_Old = sum(obs_choice *(birthNum >1))/sum(birthNum >1),
                                                           WTP = sum(log(1/(1-hospShare))),
                                                           WTP_Old = sum(log(1/(1-hospShare))*(birthNum >1)), N = .N,
                                                           NOld = sum(birthNum >1)),
                             by = c("pat_zip","year")]
        dataWTP2[, hospital := 2]
        dataSimulAgg = dataSimul[hosp_id_alt %in% c(hosp2List,"100167"), list(obs_choice = sum(obs_choice), hospShare = sum(hospShare)), by = c("maskssn","year", "birthNum", "pat_zip")]
        dataWTP12 = dataSimulAgg[, list(share = sum(obs_choice)/.N,
                                                                  share_Old = sum(obs_choice *(birthNum >1))/sum(birthNum >1),
                                                                  WTP = sum(log(1/(1-hospShare))),
                                                                  WTP_Old = sum(log(1/(1-hospShare))*(birthNum >1)), N = .N,
                                                                  NOld = sum(birthNum >1)),
                                 by = c("pat_zip","year")]
        dataWTP12[, hospital := 12]
        dataWTP3 = dataSimul[hosp_id_alt == "999999", list(share = sum(obs_choice)/.N,
                                                         share_Old = sum(obs_choice *(birthNum >1))/sum(birthNum >1),
                                                         WTP = sum(log(1/(1-hospShare))),
                                                         WTP_Old = sum(log(1/(1-hospShare))*(birthNum >1)), N = .N,
                                                         NOld = sum(birthNum >1)),
                           by = c("pat_zip","year")]
        dataWTP3[, hospital := 3]
      
        resultsList[[xj]]= rbind(dataWTP1, dataWTP2, dataWTP12, dataWTP3)
      }
  return(resultsList)
  }

doSimulationWrap = function(ITER,dataCohortInitial,dataCohortNew,demandOld,dataPSAMktList,DISTANCE,SWITCHING.COST,YEAR.START,YEAR.END) {
    dataSimulList= list(list(mean = list(),plusHalf= list(), plusOne= list(), minusHalf = list(), minusOne = list(), old = list()),
                        list(mean = list(),plusHalf= list(), plusOne= list(), minusHalf = list(), minusOne = list(), old = list()),
                        list(mean = list(),plusHalf= list(), plusOne= list(), minusHalf = list(), minusOne = list(), old = list()))
  for(xi in 1:3) {
    set.seed(123457 + ITER)
    allBirthCohorts = createInitialCohort(dataCohortInitial,dataCohortNew,demandOld,YEAR.START,YEAR.END)
    ### Need to deal w/ outside option
    dataMerged = copy(merge(allBirthCohorts,dataPSAMktList[[xi]], by = c("pat_zip"), allow.cartesian = TRUE))
    dataMerged[,chosenPrev := hosp_id_prev == hosp_id_alt]
    dataMerged[hosp_id_alt == "0",chosenPrev := FALSE]
    dataMerged[,rand := runif(1), by = c("maskssn","birthNum")]
    results = doSimulation(ITER = ITER,data = dataMerged,DISTANCE = DISTANCE[xi],SWITCHING.COST = SWITCHING.COST[xi],YEAR.START = YEAR.START,YEAR.END = YEAR.END)
    for(xj in 1:6) {
      dataSimulList[[xi]][[xj]] = results[[xj]]
    }
  }
    return(dataSimulList)
}


YEAR.START = 2015
YEAR.END = 2025

time.start = proc.time()
ITER.LIMIT = 100
simul = mclapply(1:ITER.LIMIT,
    function(x) doSimulationWrap(ITER = x,dataCohortInitial = dataCohortPSA,dataCohortNew = newMotherCohort,demandOld = demandOld,dataPSAMktList = dataPSAMktList,DISTANCE = DISTANCE, SWITCHING.COST = SWITCHING.COST,YEAR.START = YEAR.START,YEAR.END = YEAR.END),
    mc.cores = 50)
time.finish = proc.time()
time.finish - time.start

saveRDS(simul,paste0(resultsDir,"/simulListEntry.rds"))



simulResults = rbind(cbind(do.call("rbind",lapply(1:ITER.LIMIT,function(x) simul[[x]][[1]][[1]])), CF ="mean", type = "Fixed Effect"),
                     cbind(do.call("rbind",lapply(1:ITER.LIMIT,function(x) simul[[x]][[2]][[1]])), CF="mean", type = "Lagged Dep Var"),
                     cbind(do.call("rbind",lapply(1:ITER.LIMIT,function(x) simul[[x]][[3]][[1]])), CF="mean", type = "Standard Logit"),
                     cbind(do.call("rbind",lapply(1:ITER.LIMIT,function(x) simul[[x]][[1]][[2]])), CF= "plusHalf", type = "Fixed Effect"),
                     cbind(do.call("rbind",lapply(1:ITER.LIMIT,function(x) simul[[x]][[2]][[2]])), CF = "plusHalf", type = "Lagged Dep Var"),
                     cbind(do.call("rbind",lapply(1:ITER.LIMIT,function(x) simul[[x]][[3]][[2]])), CF = "plusHalf", type = "Standard Logit"),
                     cbind(do.call("rbind",lapply(1:ITER.LIMIT,function(x) simul[[x]][[1]][[3]])), CF = "plusOne", type = "Fixed Effect"),
                     cbind(do.call("rbind",lapply(1:ITER.LIMIT,function(x) simul[[x]][[2]][[3]])), CF = "plusOne", type = "Lagged Dep Var"),
                     cbind(do.call("rbind",lapply(1:ITER.LIMIT,function(x) simul[[x]][[3]][[3]])), CF = "plusOne", type = "Standard Logit"),
                     cbind(do.call("rbind",lapply(1:ITER.LIMIT,function(x) simul[[x]][[1]][[4]])), CF= "minusHalf", type = "Fixed Effect"),
                     cbind(do.call("rbind",lapply(1:ITER.LIMIT,function(x) simul[[x]][[2]][[4]])), CF = "minusHalf", type = "Lagged Dep Var"),
                     cbind(do.call("rbind",lapply(1:ITER.LIMIT,function(x) simul[[x]][[3]][[4]])), CF = "minusHalf", type = "Standard Logit"),
                     cbind(do.call("rbind",lapply(1:ITER.LIMIT,function(x) simul[[x]][[1]][[5]])), CF = "minusOne", type = "Fixed Effect"),
                     cbind(do.call("rbind",lapply(1:ITER.LIMIT,function(x) simul[[x]][[2]][[5]])), CF = "minusOne", type = "Lagged Dep Var"),
                     cbind(do.call("rbind",lapply(1:ITER.LIMIT,function(x) simul[[x]][[3]][[5]])), CF = "minusOne", type = "Standard Logit"),
                     cbind(do.call("rbind",lapply(1:ITER.LIMIT,function(x) simul[[x]][[1]][[6]])), CF = "old", type = "Fixed Effect"),
                     cbind(do.call("rbind",lapply(1:ITER.LIMIT,function(x) simul[[x]][[2]][[6]])), CF="old", type = "Lagged Dep Var"),
                     cbind(do.call("rbind",lapply(1:ITER.LIMIT,function(x) simul[[x]][[3]][[6]])), CF="old", type = "Standard Logit")
)


simulResultsAgg = simulResults[,lapply(.SD,mean), by = c("pat_zip","year","CF","type","hospital")]

simulResultsCasted = reshape(simulResultsAgg, timevar = "hospital", idvar = c("pat_zip","year","CF","type"), direction = "wide")

simulResultsNew = simulResultsCasted[, list(WTPM = sum(WTP.12), WTPM_Old = sum(WTP_Old.12)), by = c("year","CF","type")]

simulResultsOld = simulResultsCasted[CF == "old",]
simulResultsOld = simulResultsOld[, list(WTPO  = sum(WTP.1 + WTP.2), WTPO_Old = sum(WTP_Old.1 + WTP_Old.2)), by = c("year","type")]


simulResultsMerged = merge(simulResultsNew, simulResultsOld, by = c("year", "type"))


simulResultsMerged[,WTPC := (WTPM - WTPO)/WTPO]
simulResultsMerged[,WTPC_Old := (WTPM_Old - WTPO_Old)/WTPO_Old]

simulResultsMerged2015 = simulResultsMerged[year == 2015,]

simulResultsMerged2015 = simulResultsMerged2015[, list(type, CF, WTPC, WTPC_Old)]
setnames(simulResultsMerged2015, c("WTPC", "WTPC_Old"), c("WTPC_15","WTPC_Old_15"))

simulResultsMerged = merge(simulResultsMerged, simulResultsMerged2015, by = c("type", "CF"))

simulResultsMerged[, WTPC_Change := WTPC - WTPC_15]
simulResultsMerged[, WTPC_Old_Change := WTPC_Old - WTPC_Old_15]


simulResultsMerged[, CF := factor(CF, levels = c("old","minusOne","minusHalf","mean", "plusHalf", "plusOne"), 
                                  labels = c("No Entry", "Minus One", "Minus Half", "Mean", "Plus Half", "Plus One"))]

qplot(data = simulResultsMerged[ type != "Standard Logit",], x = year - 2015, y = WTPC_Change * 100, color = CF, geom = "line") +
  facet_grid(type ~ .) + theme_bw() + geom_line(size = 2) + theme_bw() +
  #    +
  theme(axis.title.x = element_text(face = "bold", size = 14, vjust = .1), axis.text.x = element_text(size = 10)) +
  theme(axis.title.y = element_text(face = "bold", size = 14, vjust = .1), axis.text.y = element_text(size = 10)) + 
  theme(panel.grid.major = element_blank(),
        panel.grid.minor = element_blank() ) + guides(color = guide_legend(nrow =2, byrow = TRUE )) +  theme(strip.text.y = element_text(size = 12, face = "bold")) +
  theme(legend.position = "bottom", legend.box = "vertical") + 
  theme(legend.title = element_blank(), legend.text = element_text(size = 12)) +
  labs(x = "Year", y = "Percent Change in WTP") + scale_x_continuous(breaks = c(0, 2, 4, 6, 8, 10)) +scale_y_continuous(breaks = pretty_breaks(n=6)) 

VAR = "WTP"
CF_VALUE = "Change"

ggsave(filename = paste0(resultsDir,"/Switch/Entry_",CF_VALUE,"_",VAR,"_Pres",".png"), width = widthPres,
       height = heightPres, dpi = dpiPlot, units = unitsPres, type = "cairo")
ggsave(filename =  paste0(resultsDir,"/Switch/Entry_",CF_VALUE,"_",VAR,"_Paper",".png"), width = widthPaper,
       height = heightPaper, dpi = dpiPlot, units = unitsPaper, type = "cairo")


qplot(data = simulResultsMerged[ type != "Standard Logit",], x = year - 2015, y = WTPC_Change * 100, color = CF, geom = "line") +
  facet_grid(type ~ .) + theme_bw() + geom_line(size = 2) + theme_bw() +
  #    +
  theme(axis.title.x = element_text(face = "bold", size = 14, vjust = .1), axis.text.x = element_text(size = 10)) +
  theme(axis.title.y = element_text(face = "bold", size = 14, vjust = .1), axis.text.y = element_text(size = 10)) + 
  theme(panel.grid.major = element_blank(),
        panel.grid.minor = element_blank() ) + guides(color = guide_legend(nrow =2, byrow = TRUE )) +  theme(strip.text.y = element_text(size = 12, face = "bold")) +
  theme(legend.position = "bottom", legend.box = "vertical") + 
  theme(legend.title = element_blank(), legend.text = element_text(size = 12)) +
  labs(x = "Year", y = "Percent Change in WTP") + scale_x_continuous(breaks = c(0, 2, 4, 6, 8, 10)) +scale_y_continuous(breaks = pretty_breaks(n=6)) + scale_color_grey()

VAR = "WTP"
CF_VALUE = "Change"

ggsave(filename = paste0(resultsDir,"/Switch/Entry_",CF_VALUE,"_",VAR,"_PresBW",".png"), width = widthPres,
       height = heightPres, dpi = dpiPlot, units = unitsPres, type = "cairo")
ggsave(filename =  paste0(resultsDir,"/Switch/Entry_",CF_VALUE,"_",VAR,"_PaperBW",".png"), width = widthPaper,
       height = heightPaper, dpi = dpiPlot, units = unitsPaper, type = "cairo")

qplot(data = simulResultsMerged[year == 2015,], x = CF, y = WTPC * 100, fill = type, geom = "bar") + 
  theme_bw() + geom_bar(stat = "identity", position = "dodge")  +
  theme(axis.title.x = element_text(face = "bold", size = 14, vjust = .1), axis.text.x = element_text(size = 10)) +
  theme(axis.title.y = element_text(face = "bold", size = 14, vjust = .1), axis.text.y = element_text(size = 10)) + 
  theme(panel.grid.major = element_blank(),
        panel.grid.minor = element_blank() ) +
  theme(legend.position = "bottom", legend.box = "vertical") +  guides(fill = guide_legend(nrow =2, byrow = TRUE )) +
  theme(legend.title = element_blank(), legend.text = element_text(size = 12)) +
  labs(x = "Counterfactual Scenario", y = "Percent Change in WTP") 


VAR = "WTP"
CF_VALUE = "Level"

ggsave(filename = paste0(resultsDir,"/Switch/Entry_",CF_VALUE,"_",VAR,"_Pres",".png"), width = widthPres,
       height = heightPres, dpi = dpiPlot, units = unitsPres, type = "cairo")
ggsave(filename =  paste0(resultsDir,"/Switch/Entry_",CF_VALUE,"_",VAR,"_Paper",".png"), width = widthPaper,
       height = heightPaper, dpi = dpiPlot, units = unitsPaper, type = "cairo")


qplot(data = simulResultsMerged[year == 2015,], x = CF, y = WTPC * 100, fill = type, geom = "bar") + 
  theme_bw() + geom_bar(stat = "identity", position = "dodge")  +
  theme(axis.title.x = element_text(face = "bold", size = 14, vjust = .1), axis.text.x = element_text(size = 10)) +
  theme(axis.title.y = element_text(face = "bold", size = 14, vjust = .1), axis.text.y = element_text(size = 10)) + 
  theme(panel.grid.major = element_blank(),
        panel.grid.minor = element_blank() ) +
  theme(legend.position = "bottom", legend.box = "vertical") +  guides(fill = guide_legend(nrow =2, byrow = TRUE )) +
  theme(legend.title = element_blank(), legend.text = element_text(size = 12)) +
  labs(x = "Counterfactual Scenario", y = "Percent Change in WTP")  + scale_fill_grey()


VAR = "WTP"
CF_VALUE = "Level"

ggsave(filename = paste0(resultsDir,"/Switch/Entry_",CF_VALUE,"_",VAR,"_PresBW",".png"), width = widthPres,
       height = heightPres, dpi = dpiPlot, units = unitsPres, type = "cairo")
ggsave(filename =  paste0(resultsDir,"/Switch/Entry_",CF_VALUE,"_",VAR,"_PaperBW",".png"), width = widthPaper,
       height = heightPaper, dpi = dpiPlot, units = unitsPaper, type = "cairo")

simulResultsEntry = simulResultsCasted[, list(shareEntry = sum(N.3 * share.3)/sum(N.3),
                                              shareEntry_Old = sum(NOld.3 * share_Old.3)/sum(NOld.3)), by = c("year","CF","type")]


simulResultsEntry2015 = simulResultsEntry[year == 2015,]

simulResultsEntry2015 = simulResultsEntry2015[, list(type, CF, shareEntry, shareEntry_Old)]
setnames(simulResultsEntry2015, c("shareEntry", "shareEntry_Old"), c("shareEntry_15","shareEntry_Old_15"))

simulResultsEntry = merge(simulResultsEntry, simulResultsEntry2015, by = c("type", "CF"))

simulResultsEntry[, shareEntry_Change := shareEntry - shareEntry_15]
simulResultsEntry[, shareEntry_Old_Change := shareEntry_Old - shareEntry_Old_15]

simulResultsEntry[, CF := factor(CF, levels = c("old","minusOne","minusHalf","mean", "plusHalf", "plusOne"), 
                                  labels = c("No Entry", "Minus One SD", "Minus Half SD", "Mean", "Plus Half SD", "Plus One SD"))]

qplot(data = simulResultsEntry[year == 2015 & CF != "No Entry",], x = CF, y = shareEntry*100, fill = type, geom = "bar")+ 
  theme_bw() + geom_bar(stat = "identity", position = "dodge")  +
  theme(axis.title.x = element_text(face = "bold", size = 14, vjust = .1), axis.text.x = element_text(size = 10)) +
  theme(axis.title.y = element_text(face = "bold", size = 14, vjust = .1), axis.text.y = element_text(size = 10)) + 
  theme(panel.grid.major = element_blank(),
        panel.grid.minor = element_blank() ) + 
  theme(legend.position = "bottom", legend.box = "vertical") + guides(fill = guide_legend(nrow =2, byrow = TRUE )) +
  theme(legend.title = element_blank(), legend.text = element_text(size = 12)) +
  labs(x = "Counterfactual Scenario", y = "Share") 


VAR = "Share"
CF_VALUE = "Level"

ggsave(filename = paste0(resultsDir,"/Switch/Entry_",CF_VALUE,"_",VAR,"_Pres",".png"), width = widthPres,
       height = heightPres, dpi = dpiPlot, units = unitsPres, type = "cairo")
ggsave(filename =  paste0(resultsDir,"/Switch/Entry_",CF_VALUE,"_",VAR,"_Paper",".png"), width = widthPaper,
       height = heightPaper, dpi = dpiPlot, units = unitsPaper, type = "cairo")


qplot(data = simulResultsEntry[year == 2015 & CF != "No Entry",], x = CF, y = shareEntry*100, fill = type, geom = "bar")+ 
  theme_bw() + geom_bar(stat = "identity", position = "dodge")  +
  theme(axis.title.x = element_text(face = "bold", size = 14, vjust = .1), axis.text.x = element_text(size = 10)) +
  theme(axis.title.y = element_text(face = "bold", size = 14, vjust = .1), axis.text.y = element_text(size = 10)) + 
  theme(panel.grid.major = element_blank(),
        panel.grid.minor = element_blank() ) + 
  theme(legend.position = "bottom", legend.box = "vertical") + guides(fill = guide_legend(nrow =2, byrow = TRUE )) +
  theme(legend.title = element_blank(), legend.text = element_text(size = 12)) +
  labs(x = "Counterfactual Scenario", y = "Share") + scale_fill_grey()


VAR = "Share"
CF_VALUE = "Level"

ggsave(filename = paste0(resultsDir,"/Switch/Entry_",CF_VALUE,"_",VAR,"_PresBW",".png"), width = widthPres,
       height = heightPres, dpi = dpiPlot, units = unitsPres, type = "cairo")
ggsave(filename =  paste0(resultsDir,"/Switch/Entry_",CF_VALUE,"_",VAR,"_PaperBW",".png"), width = widthPaper,
       height = heightPaper, dpi = dpiPlot, units = unitsPaper, type = "cairo")

qplot(data = simulResultsEntry[ type != "Standard Logit",], x = year - 2015, y =shareEntry_Change*100, color = CF, geom = "line") +
  facet_grid(type ~ .) + theme_bw() + geom_line(size = 2) + theme_bw() +
  #    +
  theme(axis.title.x = element_text(face = "bold", size = 14, vjust = .1), axis.text.x = element_text(size = 10)) +
  theme(axis.title.y = element_text(face = "bold", size = 14, vjust = .1), axis.text.y = element_text(size = 10)) + 
  theme(panel.grid.major = element_blank(),
        panel.grid.minor = element_blank() ) +
  theme(legend.position = "bottom", legend.box = "vertical") + guides(color = guide_legend(nrow =2, byrow = TRUE )) +  theme(strip.text.y = element_text(size = 12, face = "bold")) +
  theme(legend.title = element_blank(), legend.text = element_text(size = 12)) +
  labs(x = "Year", y = "Change in Share") + scale_x_continuous(breaks = c(0,2,4,6,8,10)) +scale_y_continuous(breaks = pretty_breaks(n=6)) 


VAR = "Share"
CF_VALUE = "Change"

ggsave(filename = paste0(resultsDir,"/Switch/Entry_",CF_VALUE,"_",VAR,"_Pres",".png"), width = widthPres,
       height = heightPres, dpi = dpiPlot, units = unitsPres, type = "cairo")
ggsave(filename =  paste0(resultsDir,"/Switch/Entry_",CF_VALUE,"_",VAR,"_Paper",".png"), width = widthPaper,
       height = heightPaper, dpi = dpiPlot, units = unitsPaper, type = "cairo")


qplot(data = simulResultsEntry[ type != "Standard Logit",], x = year - 2015, y =shareEntry_Change*100, color = CF, geom = "line") +
  facet_grid(type ~ .) + theme_bw() + geom_line(size = 2) + theme_bw() +
  #    +
  theme(axis.title.x = element_text(face = "bold", size = 14, vjust = .1), axis.text.x = element_text(size = 10)) +
  theme(axis.title.y = element_text(face = "bold", size = 14, vjust = .1), axis.text.y = element_text(size = 10)) + 
  theme(panel.grid.major = element_blank(),
        panel.grid.minor = element_blank() ) +
  theme(legend.position = "bottom", legend.box = "vertical") + guides(color = guide_legend(nrow =2, byrow = TRUE )) +  theme(strip.text.y = element_text(size = 12, face = "bold")) +
  theme(legend.title = element_blank(), legend.text = element_text(size = 12)) +
  labs(x = "Year", y = "Change in Share") + scale_x_continuous(breaks = c(0,2,4,6,8,10)) +scale_y_continuous(breaks = pretty_breaks(n=6))  + scale_color_grey()


VAR = "Share"
CF_VALUE = "Change"

ggsave(filename = paste0(resultsDir,"/Switch/Entry_",CF_VALUE,"_",VAR,"_PresBW",".png"), width = widthPres,
       height = heightPres, dpi = dpiPlot, units = unitsPres, type = "cairo")
ggsave(filename =  paste0(resultsDir,"/Switch/Entry_",CF_VALUE,"_",VAR,"_PaperBW",".png"), width = widthPaper,
       height = heightPaper, dpi = dpiPlot, units = unitsPaper, type = "cairo")

